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ABSTRACT 


An  algorithm  to  solve  linear  programming  problems  is  presented 
which  is  based  on  Karmarkar's  projective  method.  The  algorithm 
includes  a  practical  method  to  project  a  general  linear  programming 
problem  onto  a  unit  simplex  and  eliminates  the  a  priori  need  to  know 
the  optimal  value  of  the  objective  function.  The  implementation 
conserves  sparsity.  The  key  part  of  the  implementation  is  the  solution 
of  a  linear  least- squares  problem  to  find  an  improving  direction:  a 
direct  and  an  iterative  method  are  implemented  to  solve  this  problem. 
The  direct  method  employs  the  minimum-degree  heuristic  to  reorder  the 
system  of  normal  equations,  and  thus  conserve  sparsity  during  the 
following  Cholesky  factorization.  The  iterative  method  uses  the  incom¬ 
plete  Cholesky  factor  of  the  normal  equation  matrix  as  a  preconditioner 
for  conjugate  gradient  iterations  which  are  performed  implicitly  on  the 
preconditioned  matrix.  The  study  concludes  with  implementation 
remarks,  and  computational  results. 
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I.  INTRODUCTION 


In  recent  papers  Karmarkar  [Refs.  1,2]  has  presented  a  new 
method  for  the  solution  of  linear  programming  (LP)  problems.  His  new 
solution  technique  moves  from  some  feasible  starting  point  across  the 
interior  region  of  a  polytope  that  is  defined  by  the  problem  constraints. 
He  shows  that  the  number  of  steps  to  find  an  optimal  solution  with  his 
technique  is  polynomially  bounded. 

In  contrast,  the  simplex  algorithm  which,  is  widely  used  for  the 
solution  of  LP  problems  finds  an  optimal  solution  by  moving  from  vertex 
to  vertex  on  the  poly  tope.  This  is  known,  in  the  worst  case,  to  require 
an  exponential  number  of  steps  [Ref.  3]. 

The  polynomial  bound  of  the  projective  algorithm  makes  the  new 
solution  technique  very  appealing  to  researchers.  The  theory  of  the 
new  algorithm  seems  to  be  widely  accepted  among  experts,  while 
Karmarkar's  claim  that  his  algorithm  is  50-100  times  faster  than  the 
simplex  method  has  met  skepticism  [Ref.  4]. 

In  this  study  a  variant  of  the  projective  method  is  implemented, 
anc  some  well  known  test  problems  are  solved. 

A.  THE  BASIC  PROJECTIVE  METHOD 

Following  Karmarkar  [Ref.  1:  p.  4],  the  number  of  steps  of  the 
algorithm  depends  on  R/r,  where  R  is  the  radius  of  the  sphere 
circumscribing  the  polytope,  and  r  the  radius  of  the  inscribed  sphere. 

Assume  a  general  LP  of  the  form 

T 

Min  c  x 

s.t.  A  x  =  b  (LP1) 

x  >  0  . 


With  the  assumption  that  the  sum  of  its  variables  has  an  upper  bound, 
and  with  the  proper  scaling  of  variables,  a  convexity  constraint 


can  be  added  to  LP1.  This  transformation  maps  the  LP  onto  a  unit 
simplex  S  whose  center  is  at  aQ  =  (1/n,  .  .  . ,  1/n) .  Then,  a  transforma¬ 
tion  is  performed  such  that  a  feasible  interior  starting  point  is  mapped 
onto  aQ. 

If  B(aQ,r)  is  the  largest  sphere  with  center  aQ  that  can  be 
inscribed  into  the  simplex  S,  and  B(aQ,R)  is  the  smallest  sphere 
circumscribing  S,  then 

R/r  =  n.  (1.2) 

By  restricting  a  solution  x  to  remain  inside  the  largest  inscribed  sphere 
B(aQ,r),  the  method  achieves  in  one  iteration  a  reduction  in  the  differ¬ 
ence  between  the  current  objective  value  and  the  optimal  objective  value 
by  a  factor  of  (1  -  1/n).  Following  Shanno  [Ref.  5J  a  simple  proof  is: 


Let 

f'  =  min  c^x,  x  €  S,  Ax  =  b,  (1.3) 

X  =  min  cTx,  x  €  B(a0,r),  Ax  =  b,  lTx  =1,  (1.4) 

T  =  min  cTx,  x  e  B (aQ, R) ,  Ax  =  b,  lTx  =  1.  (1.5) 

Then, 

cTa0  -  i  <  cTa0  -  f*  <  cTa0  -  f ,  (1.6) 

and  with  equation  1.2 

cTa0  -  f  =  n(cTa0  -  _f ) .  (1.7) 

From  that 

(L  -  0/(cTa„  -  f*)  <  (1  -  1/n). 


(1.8) 


If  a  linear  objective  function  could  be  maintained  at  each  iteration, 
it  follows  that  an  upper  bound  on  the  number  of  steps  required  to  find 
an  optimal  solution  is  of  O(n).  Unfortunately,  the  projective  transfor¬ 
mations  needed  to  continue  the  algorithm  result  in  a  nonlinear  objective 
function. 

B.  DERIVATION  OF  THE  ALGORITHM 
1.  The  Canonical  Form 

Suppose  we  have  a  linear  programming  problem  of  the  form 
LP1.  Karmarkar's  method  requires  that  this  LP  be  transformed  into  the 
following  canonical  form, 

T 

Min  c  x 

s.t.  A  x  =  0  (LP2) 

lTx  =  1 

x  >  0, 

where  the  optimal  solution  value  is  zero.  With  the  assumption  that  the 

sum  of  the  variables  of  an  LP  is  bounded  above  and  subsequent  scaling 
by  this  bound,  a  convexity  constraint  (equation  1.1)  can  be  added  to 
LP1.  In  practice  this  can  be  a  problem  because  choosing  too  large  an 
upper  bound  may  cause  numerical  problems. 

LP2  requires  that  the  nonhomogeneous  system  of  equations 
Ax  =  b  be  made  homogeneous.  Karmarkar  [Ref.  1:  p.34]  proposes  a 
transformation  that  would  transform  the  i-th  equation  ajTx  =  bj  to 

paij  -  bi)*j  =  °'  (1-9) 

This  transformation  has  the  disadvantage  that  when  b  is  dense  the 
sparsity  of  A  will  be  lost.  Karmarkar  requires  the  optimal  solution 
value  to  LP2  to  be  zero  together  with  a  special  stopping  criterion 
(equation  1.26),  to  prove  the  polynomial  bound.  To  achieve  the  zero 
objective  value  the  optimal  objective  function  value  f  of  the 


untransformed  LP  has  to  be  known  in  advance,  and  the  following  trans¬ 
formation  made 


cTx  -  f"  =  cTx  -  f"lTx  =  (c  -fl)Tx  . 


(1.10) 


Karmarkar  concludes  [Ref.  2:  p.  387]  that  if  the  minimal  objective  value 
determined  by  the  projective  algorithm  is  not  equal  to  zero,  the  original 
problem  must  be  either  infeasible  or  unbounded.  Another  method  for 
making  the  transformation  from  LP1  to  LP2  is  discussed  in 
Chapter  II.  B. 

2 .  The  Projective  Transformation 

Let  x°  >  0  be  a  known  feasible  solution  to  LP2.  The  invertible 
projective  transformation 


y  = 


D"1x 


1TD'1x 


(1.11) 


T 

maps  any  x,  such  that  lAx  =  1  and  x  >  0,  onto  y  such  that 

T 

1  y  =  1  and  y  >  0.  D  is  a  diagonal  matrix  whose  entries  are 

<xl° . xn°)  • 

The  transformation  maps  the  LP2  unit  simplex  in  x- space  onto 

another  unit  simplex  in  y-space.  The  point  x°  is  mapped  onto 

o  T 

y  =  l/n  1  ,  the  center  of  the  unit  simplex  in  y-space.  The  inverse  of 
the  transformation  (equation  1.11)  is  given  by 


x  = 


D  y 

1TD  y 


(1.  12) 


After  the  transformation,  LP2  can  now  be  restated  as 


cTP  y 
1TD  y 


s.t.  A  D  y  =  0 


(LP3) 


iT  y 


>  o  . 


Strictly  speaking  LP3  is  not  a  linear  program  because  its  objective 
function  is  a  rational  function  of  y.  Karmarkar  [Ref.  1:  page  18ff] 
shows  that  it  is  sufficient  to  consider  a  linear  approximation  to  the 
objective  function  of  LP3,  which  leads  to 


Min  c 1 D  y 
s.t.  A  D  y  =  0 

1T  y  =1 


(LP4) 


3.  Optimization  Over  a  Sphere 


A  solution  to  LP4  is  now  restricted  to  lie  within  a  sphere  with 


center  at  y  and  radius  «r,  where 


r  =  1/ (n(n- 1) ) ' 1/2  . 


(1.13) 


r  is  the  radius  of  the  largest  sphere  that  can  be  inscribed  into  the  unit 
simplex,  and  a  is  a  constant  such  that  0<a<l.  a  provides  a  margin 
which  ensures  that  the  algorithm  doesn't  select  a  point  outside  the 
sphere  due  to  round-off  error.  By  convexity,  an  optimal  solution  will 
occur  at  the  boundary  of  the  sphere.  Thus,  the  additional  constraint 
can  be  stated  as  an  equality.  Now  we  have,  the  following  version  of  the 
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VV-* 


T 

Min  c  D  y 
s.t.  A  D  y 


(y-y°)T(y-y°) 


(LP5) 


=  0C2r2 


Before  the  problem  can  be  solved,  one  more  transformation 
that  moves  the  center  of  the  sphere  to  the  origin  is  useful.  Let 


y  =  y  +  y°,  and  eliminate  the  constant  terms  ADy°  =  0,  l^y0  =  1  and 
c^Dy0.  For  convenience,  define  B  by 


(1.14) 


and  the  gradient  c  of  the  objective  function  of  LP5  by 


c  =  cTD  . 


(1.15) 


Then  LP5  can  be  restated  as. 


Min  cTy 


(LP6) 


=  *2r2 


In  order  to  solve  LP6,  note  that  an  initial  feasible  solution  to 
LP6  is  y  =  0  (which  corresponds  to  y  =  y°  in  LP5).  This  is  also  the 

_ rp _  o  2 

center  of  the  sphere  y  y  =  a  r.  An  optimal  solution  to  LP6  can  be 
obtained  by  finding  a  direction  of  maximum  rate  of  ascent  c  that  is 
feasible  with  respect  to  By  =  0,  and  moving  in  direction  -c  (maximum 
rate  of  descent)  a  distance  «r  from  y  =  0  to  the  boundary  of  the 
sphere . 

A  feasible  direction  of  maximum  ascent  is  found  by  orthogo¬ 
nally  projecting  c  onto  the  null  space  of  B  (see  [Ref.  1:  page  17]  ), 
i.e.,  the  following  problem  has  to  be  solved  in  terms  of  cD 
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Mm  cTc  -  2cTCp  +  cpTcp  =  Min  ( ||  5  -  cp||  2)2  (1.16) 


s.t.  BCp  =  0  . 


Define  the  Lagrangian  L(cD,  X)  by 


L(Cp,X)  =  cTc  -  2cTcp  +  CpTcp  -  xTBcp. 


(1.17) 


First-order  optimality  conditions  for  the  Lagrangian  yield 


-2c  +  2c  -  BiX  =  0, 


(1.18) 


Bcp  =  0. 


(1-19) 


After  multiplying  by  B  and  dropping  2Bcp=0  (equation  1.19)  we  get 


-2Bc  =  BBT\  . 


(1.20) 


Assuming  (BBT)~1  exists,  we  can  solve  for  A 


A  =  -2(BBT)'1Bc. 


(1.21) 


Substituting  equation  1.21  into  equation  1.18  gives 


c  =  c  -  BT(BBT)'1Bc. 


(1.22) 


The  direction  of  maximum  rate  of  ascent  is  then 


c  =  V1  cPii  • 


(1.2 3) 


Thus,  the  optimal  solutions  to  LP6  and  LP5  are 


y  =  -  <*rc, 


(1.24) 
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iM  v_r  Krw’wr*?'  «.t."  v  »_ 


rr*;v 


'1 


and 

y1  =  y°  -  A rc  (1-25) 

respectively. 

Finding  Cp  is  the  key  part  of  Karmarkar's  method  because  it 
involves  the  major  portion  of  the  computational  work.  Solving  equation 
1.22  for  Cp  can  be  viewed  as  solving  a  linear  least- squares  problem 
(see  Chapter  II.C.). 

With  the  optimal  solution  to  LP5  in  y-space,  the  method 
proceeds  by  transforming  that  solution  back  into  x-space  by  use  of 
equation  1.12.  Next  it  defines  a  new  matrix  D  and  iterates  until  it 
reaches  the  stopping  criterion  [Ref.  1:  p.  14] 

(cTxk)/(cTx°)  <  2_c*  (1.26) 

rn 

where  q  is  a  termination  parameter.  Note  that  c  x  =  0  at  optimality. 
Table  1  shows  an  outline  of  Karmarkar's  basic  method. 
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TABLE  1 
ALGORITHM  1 

Problem  Size . n 

Coefficient  Matrix . '.A 

n 

Cost  Function . cJ 

Initial  Feasible  Solution. . .xl 

Termination  Parameter . q 

Feasibility  Parameter . a 


i  k  =  1 

x  =  x° 

f°  =  CTX° 
f  =  oo 

r  =  l/(n(n-l) ) '1/2 

While  (f/f°  >  2' 3) 

D  =  diag  (x) 
y  =  (D'1x)/(1TD'1x) 

c  =  c^D 

F  A  D  1 


cp  =  c  -  BT(BBT)'1Bc 

6  =  cp/HcpH 

y  =  y  -  a  rc 
x  =  ( Dy )/( lTDy ) 

f  =  cTx 
k  =  k  +  1 

End (While) 

t  Solution . x 

Objective  Function  Value... f 
Iteration  Count . k 


II.  A  VARIANT  QF  THE  PROJECTIVE  METHOD 


A.  INTRODUCTION 

The  following  sections  outline  a  variant  of  the  projective  method 
that  was  proposed  by  Wood  [Ref.  6].  It  has  a  practical  method  of 
bringing  an  LP  into  the  canonical  form,  and  uses  the  non-linear  objec¬ 
tive  function  of  LP3  to  find  a  gradient  c'  which  corresponds  to  c  of 
equation  1.16.  Similar  to  the  simplex  method,  the  proposed  algorithm 
uses  a  ratio  test  to  determine  a  feasible  step  length. 

Other  practical  features  of  the  proposed  method  include  the  relax¬ 
ation  of  the  requirement  to  know  the  optimal  objective  value  in  advance, 
and  exploitation  of  sparsity  of  A.  The  implementation  is  especially 
concerned  with  controlling  fill-in  during  the  solution  of  the  normal 
equations. 

B.  DERIVATION  OF  THE  MODIFIED  ALGORITHM 

Given  an  LP  problem  of  the  form  LP1,  we  use  a  single  artificial 
variable  xn  +  1  to  attain  initial  feasibility: 

Min  (cT,M)(x,xn  +  1) 

s.t.  Ax-(Ax°-b)xn  +  1  =  b  (LP7) 

x,xn+1  >  0  . 

where  M  is  the  cost  for  the  artificial  variable  and  x°  is  the  initial  solu¬ 
tion  with 

(x°,x0n+i)  =  1T  .  (2.1) 

The  transformation  has  added  one  more  column  to  the  problem. 

To  get  a  homogeneous  right-hand  side  in  LP7  we  introduce  an 
additional  variable,  and  apply  the  following  projective  transformation, 
whose  inverse  is  given  by 
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V  V  V 

I  *  *  »  * 


'V. 


>;-r '  w\  TT  w  ■;  f  ‘,  » 


t  r  v”  ■- 


r  v *7 


^  wj-»  »j-  • 


MV 


a". 


(x’xn+l)  =  (n+2)(x',x'n+1)/x'n+2. 


(2.2) 


(x1 ,  x'n+ 1 )  =  (n+2)(x,xn+1)/(lT(x,xn  +  1)+l)  (2.3) 

x'n+2  =  (n+2)/ (lT(x, xR+1)  +  1)  .  (2.4) 

Let  the  artificial  column  be  denoted  by  a,  i.e.  a  =  -(Ax°-b).  The 
transformed  problem  can  now  be  restated,  with  the  exception  of  the 
objective  function,  in  Karmarkar's  canonical  form, 


Min 

(c,M,0)  (x',x'n  +  1  ,x'n+2) 

x'n+2 

s.t. 

(A, a,  -b)(x',x'n+1,x'n+2) 

=  0 

lT(x',x'n  +  1,x'n+2) 

=  n+2 

(x',x'n+1,x'n+2) 

>  0. 

(LP8) 


The  above  transformation  adds  yet  another  column  to  the  problem,  but 
has  the  advantage  that  it  doesn't  change  the  sparsity  of  A,  as  opposed 
to  Karmarkar's  proposal.  Rather  than  projecting  the  LP  problem  onto  a 
unit  simplex,  it  is  projected  onto  an  (n+2) -simplex.  This  is  done  to 
improve  numerical  stability. 

The  projective  transformation,  equation  1.12,  is  applied  to  LP8 
which  gives 


Min  (cT,M,0)D(y,yn  +  1,yn+2)/dn+2yn+2 

s.t.  (A,a,-b)D(y,yn  +  1,yn+2)  =  0  (LP9) 

iT(y,yn+1,yn+2)  =  n+2 

(y,yn+1,yn+2)  ^  o  . 

The  gradient  (compare  with  equation  1.16)  of  the  objective  function  of 
LP9,  evaluated  at  the  initial  feasible  starting  point  y°  =  1^  is  propor¬ 
tional  to 
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-C-'-X)  . 


(2.5) 


c'  =  (c1d1,...,cndn,Mdn+1, 

The  step  of  normalizing  Cp  in  Algorithm  1  is  replaced  by  a  ratio 
test.  Let 

jmax  =  argmax{(Cp)j}.  (2.6) 

This  allows  the  algorithm  to  make  a  step  outside  the  inscribed  sphere, 
but  maintains  feasibility  by  restricting  a  solution  to  lie  inside  the 
simplex.  Then,  the  update  of  y  becomes 

y  -  1  "  ^cp/(cp)jmax  »  (2-7) 

where  p  is  a  parameter  to  maintain  feasibility  such  that  0<p<l.  Table  2 
shows  the  modified  algorithm. 

C.  THE  LINEAR  LEAST- SQUARES  PROBLEM 

Computation  of  the  projected  gradient  Cp  during  every  iteration 
accounts  for  most  of  the  computational  workload  in  any  algorithm  based 
on  Karmarkar's  method.  Solving  equation  1.22  can  be  viewed  as  solving 
the  following  linear  least-squares  problem, 

Min  (  ||c'  -  BTAII2)2  (2.8) 

T 

where  B  is  an  (n  +  2)xm  matrix  with  m<(n+2)  assumed. 

T 

If  rank(B  )=m,  then  the  solution  to  (2.8)  is  given  by  the  solution 
to  the  system  of  normal  equations 

BBTA  =  Be’  .  (2.9) 

The  projected  gradient  Cp  is  then  the  residual  vector  of  the  least- 
squares  problem  (2.8),  i.e.. 


TABLE  2 
ALGORITHM  2 


Problem  Size . n 

Coefficient  Matrix . A 

T 

Cost  Function . c 

Initial  Feasible  Solution. . .x° 

Termination  Parameter . t 

Feasibility  Parameter . P 


k  = 
x  = 
f°  = 

While  (  f  °  - 


End (While ) 


1 

x° 

oo 


f(x)  >  t  ) 
D 
f° 
y 

c' 

B 


diag  (x) 
f  (x) 

( n+2 )  (D'^xJ/u'1^ 

^yf(y) 

r  a  d  i 


c’  -  BT(BBT)  ^Bc’ 


Wx  =  argmax(cp)j 

y  -  1  -  ? Cp/ ( Op )  -jmax 

x  =  (n+2)(Dy)/(l1’Dy) 

k  =  k  +  1 


t  Solution . x 

Objective  Function  Value... f 
Iteration  Count . k 


BBT  is  symmetric  and  positive  definite,  given  that  BT  has  full  column 
rank. 

As  Heath  [Ref.  7:  p.  499]  points  out,  the  ideal  choice  for  solving 
the  normal  equations,  given  full  rank,  is  Cholesky  factorization.  If  BT 
is  not  of  full  column  rank  the  cross-product  matrix  will  be  singular  and 

rP  _  1  rTt  rP 

(BB  )  ceases  to  exist.  BBi  could  become  nearly  singular  if  Bx  is 
near  rank  degenerate.  Shanno  [Ref.  5:  p.  25]  shows  that  this  will 
happen  if  the  optimal  solution  is  degenerate,  i.e.  as  the  optimum  is 
approached  numerical  problems  arise  and  Cholesky  factorization  is  likely 
to  fail. 

One  problem  that  cannot  be  avoided  when  solving  the  normal  equa¬ 
tions  is  the  fact  that  the  P-condition  number  of  BB^  is  the  square  of 
that  of  B^  [Ref.  8:  p.  223],  so  that  when  B^  is  already  ill-conditioned 
it  may  be  impossible  to  find  an  accurate  solution  to  equation  2.9. 

Another  important  consideration  when  computing  BB^  is  that  the 

T  p 

sparsity  in  B  will  not  automatically  guarantee  sparsity  in  BBX.  In 
fact,  the  addition  of  variables  in  LP7  and  LP8  has  added  two  possibly 
dense  bottom  rows  into  B  .  Thus,  BB  will  be  completely  dense. 
However,  one  can  cope  with  that  by  initially  omitting  the  dense  rows  in 

rp  pi 

B  from  the  computation  of  BB  ,  and  later  updating  the  solution  to 
equation  2.9  using  procedures  similar  to  the  ones  described  in 
[Ref.  9:  p.  58-65], 

D.  IMPLEMENTATION 

Algorithm  2  has  been  implemented  in  FORTRAN  H  (Extended) 
Opt(2)  on  an  IBM  3033  AP  under  VM/CMS.  All  floating  point  arithmetic 
is  performed  in  double  precision.  The  program  is  designed  to  accept 
different  solution  modules  from  available  software  packages  for  solving 
the  linear  least- squares  problem. 

Input  data  sets  are  in  standard  MPS  format.  The  numerical  values 
of  the  non-zero  elements  of  the  constraint  matrix  A  are  stored  column¬ 
wise  in  a  real  array.  For  versatility  a  full  set  of  pointers  are  defined: 

1.  IC  column  index 

2.  R  row  index 
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3.  AP  location  of  first  non-zero  element  in  a  column 

4.  RP  location  of  last  non-zero  element  in  a  row 

5.  LINK  location  of  next  non-zero  element  (backwards) 

in  a  row. 

Brameller,  Allan  and  Hamam  [Ref.  10:  p.  104-110]  give  a  comprehensive 
discussion  of  sparse  storage  schemes. 

The  transformation  of  an  LP  problem  into  the  canonical  form  adds 
one  dense  row  and  two  dense  columns  to  the  A  matrix.  The  dense  row 
originates  from  the  convexity  constraint  equation  1.1,  the  first  dense 
column  stems  from  the  single  artificial  variable  that  is  needed  to  attain 
initial  feasibility,  and  the  second  dense  column  is  added  to  make  the 
system  of  equations  homogeneous  (see  LP8) .  The  artificial  column  is 
updated  with  the  residual  of  the  current  solution  as  long  as  the  total 
infeasibility  is  above  a  specified  threshold. 

T  T 

As  mentioned  earlier,  dense  rows  in  B  yield  BB  completely 
dense.  With  the  following  method,  this  problem  can  be  alleviated.  The 
method  applies  to  any  number  of  dense  rows  in  B  ,  but  in  this  study 
we  are  only  concerned  about  two  dense  rows,  namely  the  ones  that 
result  from  the  transformations  that  are  performed  to  get  from  LP2,  to 
LP8  via  LP7.  Consider  the  projection  problem  of  equation  1.16  in  the 


following  form, 

Min  (  ||c'  -  cp  ||2)2  (2.11) 

s.t.  BCp  =  0 

Replace  Cp  by  z,  and  let 

B  =  (Bj,  B2)  (2.12) 

c'  =  (c'j,^)  (2.13) 

Z  =  (z1,Z2)T  -  (2.14) 
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where  is  an  (mxn^),  B2  an  (mxn2)  matrix,  and  n2  is  the  number  of 
dense  columns  in  B.  Equation  2.11  now  becomes 


Min  (  ||c'2  -  z2ll  2)2  +  (H  c'i  '  ziH  2)2  (2.15) 

s.t.  B1z1  =  -B2z2 

and  when  separated,  2.15  becomes 

f  (||  c'2  -  z2 1|  2)2  +  Min  (||  c\  -  zJI  2)2  (2-16) 

Min  J  z  1  1  Zl  1  1  z 

22  l  s.t.  BjZj  =  -B2z2  . 

Solving  the  inner  minimization  first,  the  Kuhn-Tucker  conditions  yield 


c'l  '  Z1  =  B1TA  -  (2  • l7) 

and 

BjZj  =  -B2z2  .  (2.18) 

Multiplying  2.17  with  Bj  we  get, 

Blc'l  *  Blzl  =  BlBlTA  •  (2-19) 

Substituting  2.18  into  2.19  gives 

Blc'l  +  B2Z2  =  BlBlTA.  (2.20) 

Let  A0  solve 

Blc'l  =  B1B1TA  (2.21) 

and  let  Aj  solve 

(B2)j  =  ,  j  =  1 . n2  (2.22) 
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where  (B2)j  is  the  j-th  column  of  B2,  then  2.22  corresponds  to  solving 
the  matrix  equation  system 


b2  = 


(2.23) 


where  \  is  an  (mxnj)  matrix.  Then,  the  solution  to  A  for  the  inner 
minimization  is. 


A  —  An  +  X  2ir\ 


(2.24) 


Substituting  2.24  into  2.17  gives 


j'l  -  Zj  =  B^Aq  +  B1TAz2  . 


(2.25) 


To  simplify  notation,  let  h  =  B^Aq  and  H  =  B^A 
Then  equation  2 . 25  becomes 


Zj  =  c^  -  h  -  Hz2  . 


(2.26) 


Substituting  2.2G  into  2.15,  the  overall  minimization  then  becomes 


Min  (II  c'2  -  z2H  2)2  ♦  (||  h  ♦  Hz, ||  ,)2 


(2.27) 


which  is  a  simple  least-squares  problem.  The  first-order  Kuhn-Tucker 
optimality  conditions  yield, 


c'2  -  hrH  =  (I  +  HrH)z2 


(2.28) 


which  gives,  solving  for  z2, 


z2  =  (I  +  HTH)  ’ 1  (c'2  -  hTH) 


(2.29) 


With  the  solutions  to  2.29  and  2.2G  we  have  the  desired  result, 


Cp  =  z  =  (Zl,z2) 


(2.30) 


In  practice  the  following  procedure  is  followed: 


T 

1.  compute  B^B-^ 

T 

2.  factor  ,  e.g.  using  Cholesky  factorization 

3.  solve  B^B-^A  =  B^c'^,  the  solution  is  AQ 

4.  compute  h  =  B-]_T  A  Q 

5.  solve  B^Bj_  A  =  (B2)]./  the  solution  is 

6.  solve  B^B^  A  =  (62)2/  the  solution  is  A2 

T 

7.  compute  H  =  B^  A 

8.  compute  (I  +  H^H)"^-,  (note  that  H^H  is  a  2x2 
matrix  if  there  are  2  dense  rows  in  B  ) 

9.  compute  (c'2  -  h^H) 

10.  compute  Z2 

11.  compute 

T 

The  given  procedure  is  efficient  since  the  factorization  of  BB  is 
computed  only  once  and  the  same  system  is  solved  three  times  using 
different  right  hand- sides  each  time. 

As  a  stopping  criterion  for  the  algorithm  the  following  rule  is 
used, 

IF  argmax  (|x=k  -  x,k_1|)  <  t  STOP  ,  (2.31) 

j  J  J 

where  t  is  a  real  constant.  The  convergence  criterion 

(||  cp  ||/|cTx°|)  <  t  (2.32) 

mentioned  by  Lustig  [Ref.  4:  p.  12]  can  also  be  used  to  terminate  the 
algorithm . 


III.  SOLVING  THE  LINEAR  LEAST- SQUARES  PROBLEM 


A.  INTRODUCTION 

In  this  study,  the  linear  least-squares  problem  is  solved  by 
explicitly  computing  and  solving  the  normal  equations.  Although  the 
normal  equation  approach  experiences  problems  when  applied  to  ill- 
conditioned  or  near  rank  deficient  matrices,  it  behaves  acceptably  with 
sparse  and  well -conditioned  matrices  [Ref.  11]. 

The  numerical  methods  for  solving  normal  equations  generally  fall 
into  two  classes,  direct  and  iterative  methods.  One  representative  of 
each  class  is  considered  in  this  study.  Little  can  be  said  as  to  which 
class  of  methods  is  better,  except  that  direct  methods  are  more  attrac¬ 
tive  in  terms  of  computational  work,  and  iterative  methods  may  require 
less  storage  [Ref.  12:  p.  11]. 

B.  CHOLESKY  FACTORIZATION 

The  method  implemented  is  given  in  George  and  Liu  [Ref.  12],  and 
uses  Cholesky  factorization  with  a  minimum -degree  ordering  to  solve  a 
large  sparse  positive  definite  system  of  equations.  Since  BB  is  o’-.iy 
guaranteed  to  be  positive  semidefinite,  a  modification  to  the  Cholesky 
factorization  algorithm  is  considered  in  Chapter  IV. A. 3.  to  accommodate 
the  semidefinite  case. 

The  minimum- degree  algorithm  is  a  reordering  heuristic  which 
attempts  to  reduce  fill-in  during  the  factorization  phase.  The  reordering 
phase  is  entirely  symbolic;  it  amounts  to  a  symmetric  row  and  column 
permutation  of  BB  which  corresponds  to  reordering  the  columns  in 
B  .  During  this  phase,  BB  doesn't  have  to  be  computed  numerically; 
only  its  structure  has  to  determined.  Also,  the  factorization  is  first 
performed  symbolically,  thus  allowing  a  static  data  structure  for  the 
Cholesky  factor  L.  An  outline  of  the  phases  of  the  algorithm  is  given 
in  Table  3.  See  also  Heath  [Ref.  7:  p.  499]. 


TABLE  3 

MINIMUM- DEGREE  ORDERING  /  CHOLESKY  FACTORIZATION 

ALGORITHM  (MDOC) 

T 

1.  Determine  the  nonzero  structure  of  BB  . 

T  T 

2.  Find  a  permutation  matrix  P  such  that  PBBXP 
has  a  sparse  lower  triangular  Cholesky  factor 
L. 

3.  Factor  PBBTP^  symbolically  and  set  up  the  data 
structure  for  L. 

4.  Compute  PBB^p'1'  numerically. 

5.  Factor  PBB^P^  =  LL^  numerically. 

6.  Solve  Lz  =  PBb'  (back  substitution). 

7.  Solve  L  y  =  z  (forward  substitution). 


The  Minimum-Degree  Ordering  Heuristic 


The  heuristic  finds  an  ordering  of  a  symmetric  matrix  such 

that  fill-in  is  low  when  the  matrix  is  being  factored.  The  basic  idea  is, 

T 

at  each  (simulated)  factorization  step,  to  permute  the  part  of  BB 
remaining  to  be  factored  so  that  a  column  with  the  fewest  nonzeros  is  in 
the  pivot  position.  The  implementation  consists  of  six  subroutines  that 
are  given  in  George  and  Liu  [Ref.  12:  pp.  124-137].  The  subroutines 
accept  as  input  the  adjacency  graph  associated  with  BB  represented 

by  an  adjacency  structure,  and  xeturn  as  output  a  symmetric  permuta- 

T  T 

tion  of  BB  given  as  a  permutation  vector  for  the  columns  of  B  . 

T* 

Let  G=(X,E)  be  the  adjacency  graph  of  BB  ,  where  X  is 
the  set  of  nodes,  and  E  is  the  set  of  edges.  Then,  the  nodes  corre¬ 
spond  to  the  variables  of  the  least-squares  problem,  i.e.  the  columns  of 
T 

B  .  Two  nodes  x  and  y  are  said  to  be  adjacent  if  (x,y)  is  an  edge  in 
E.  The  adjacent  set  of  Y,  Y<=X  is  defined  and  denoted  by 


Adj(Y)  =  (x  £  X- Y  |  (x, y } e  E  for  some  yeY). 


(3.1) 


An  adjacency  list  for  xg  X  is  a  list  of  all  nodes  in  Adj({x}). 
Finally,  an  adjacency  structure  for  the  graph  G  is  the  set  of  adjacency 
lists  for  all  XfX  [Ref.  12:  pp.  37-41].  The  particular  adjacency 
structure  used  in  the  ordering  heuristic  stores  elements  in  each  adja¬ 
cency  list  in  contiguous  locations.  An  entry  point  array  to  the  first 
element  in  each  list  allows  access  to  the  list. 

The  ordering  heuristic  is  based  on  graph  theory,  and  involves 
the  notion  of  elimination  graphs,  quotient  graphs,  reachable  sets  and 
indistinguishable  nodes.  George  and  Liu  [Ref.  12:  pp.  92-124]  may  be 
consulted  for  more  details . 

2.  Factorization  and  Solution 

The  components  of  the  lower  triangular  Cholesky  factor  L  of 

rp 

BBa  are  computed  using  the  so  called  "inner  product  form"  algorithm 
[Ref.  12:  p.  20].  The  elements  of  L  are  given  by, 

■jj =  <eii  \!  'V)l/2  t°ri=1-2 . m  <3-2> 

Ijj  =  (ejj  -  2^  1ik1jk^/1jj  for  i=j  +  1-J+2>  •  •  •  >m  (3-3) 

T 

where  the  e^  are  the  elements  of  BB  . 

After  the  factorization  has  been  computed,  the  following  two 
linear  systems  have  to  be  solved  (see  also  Table  3): 

Lz  =  Pb'  ,  (3.4) 

and 

LTy  =  z  .  '  (3.5) 
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Solving  system  3.4  by  back  substitution  involves  the  use  of  "inner 
products"  [Ref.  12:  p.  25],  defined  by 

Zj  =  (b'j  -  l2  likxk)/lu  for  i=l,2,...,m  ,  (3.G) 

K  *  1 

where  bk  stands  for  the  i-th  element  of  the  right-hand  side  of  (3.4). 
For  this  case  L  must  be  accessed  row-by-row. 

System  3.5  is  solved  by  forward  substitution  using  "outer 
products"  [Ref.  12:  p.  26],  defined  by 

yi  =  V^i  for  i=l, 2,  .  .  .  ,m  (3.7) 

(zi+l’ •  •  ’  zm)  (zi+l'--»zm)  "  yi^i+l,i>  •  • ’^m,i^  • 

T 

For  the  latter  case,  L  must  be  accessed  row-by-row,  or  L  column-by- 
column  instead. 

C.  INCOMPLETE  CHOLESKY  FACTORIZATION 

This  method  is  an  implementation  by  Ajiz  and  Jennings  [Ref.  13]  of 
the  incomplete  Cholesky  conjugate  gradient  algorithm  (ICCG),  whose 
theory  is  given  in  [Refs.  14,15].  The  algorithm  requires  that  the 
coefficient  matrix  of  a  set  of  simultaneous  linear  equations  be  symmetric 
and  positive  definite.  It  consists  of  two  distinct  parts,  one  being  the 
Cholesky  factorization,  which  can  be  complete  or  incomplete,  and  the 
other  a  conjugate  gradient  iteration  to  solve  a  preconditioned  linear 
system,  where  the  Cholesky  factor  serves  as  the  preconditioner. 

Golub  and  Van  Loan  [Ref.  16:  pp.  373-377]  point  out  that  precon¬ 
ditioning  is  essential  for  obtaining  good  convergence  rates  with  conju¬ 
gate  gradient  methods.  The  convergence  rate  is  closely  linked  to  the 
P-condition  number,  which  is  the  ratio  of  the  largest  to  the  smallest 

T' 

eigenvalue.  It  was  mentioned  earlier  that  the  condition  number  of  BB 

T 

will  be  the  square  of  the  one  of  B  .  Ill-conditioned  problems  have 
large  condition  numbers,  and  hence  slow  convergence.  Preconditioning 
is  a  process  of  transforming  a  linear  system  so  that  its  P-condition 
number  is  improved  [Ref.  17:  p.  979]. 


Consider  again  the  original  linear  least-squares  problem  (2.9)  with 
a  generic  right-hand  side  b'  and  x  the  unknowns, 

BBTx  =  b'.  (3.8) 

Then,  equation  3.8  can  be  preconditioned  with  a  transformation  matrix 
L  giving 

L"1BBTL"Tx  =  L"1bl,  (3.9) 

or 

L~1BBTL'Ty  =  b,  (3.10) 

*T*  _  „  1 

where  y  =  Llx  and  b  =  L  1b'.  According  to  Ajiz  and  Jennings 
[Ref.  13:  p.  950]  the  ideal  choice  of  transformation  matrix  L  is  the 

T* 

Cholesky  factor  of  BB  ,  since 

L'1BBTL'T  =  I  ,  (3.11) 

pi-ovided  one  could  perform  exact  arithmetic. 

The  objectives  of  the  incomplete  factorization  phase  of  the  algo- 

T* 

rithm  are  to  transform  BB  as  close  as  possible  to  I,  and  to  reduce 
fill-in  in  the  factor  L.  This  is  accomplished  by  discarding  some  off- 
diagonal  coefficients  during  the  factorization,  whose  magnitudes  fall 
below  a  preset  threshold  limit.  The  result  of  this  operation  is  an  incom¬ 
plete  Cholesky  factor  L  that  must  satisfy 

BBT  =  LLT  :  C  ,  (3.12) 

where  C  is  the  matrix  of  elements  omitted  from  the  factorization. 
Unfortunately,  omission  of  elements  from  the  factorization  process  can 
destroy  the  positive  definiteness  property  and  hence  lead  to  a  break¬ 
down  of  the  process.  Ajiz  and  Jennings  [Ref.  13:  pp.  950-951]  have 


shown  that  introducing  diagonal  modifications  in  C  will  retain  the  posi¬ 
tive  definiteness  property.  The  matrix  C  will  be  symmetric  and  will 
have  diagonal  elements  that  are  greater  than  or  equal  to  zero.  Thus,  C 

is  a  positive  semidefinite  matrix.  Assuming  BB  to  be  positive  definite, 

T* 

adding  C  will  result  in  LL  also  being  positive  definite. 

Convergence  of  conjugate  gradient  iterations  is  only  guaranteed 
for  the  positive  definite  case.  Thus,  a  modification  to  adapt  the 
Cholesky  factorization  to  the  positive  semidefinite  case  as  with  the 
direct  method  may  cause  slow  convergence.  The  modification  is  consid¬ 
ered  in  Section  IV.  A.  3.  Table  4  gives  an  outline  of  the  ICCG 
algorithm. 


TABLE  4 

INCOMPLETE 

CHOLESKY  -  CONJUGATE  GRADIENT  (ICCG) 
ALGORITHM 

1. 

Obtain 

L,  an  incomplete  Cholesky 

factor  of  BB^. 

2  . 

Solve 

Lb  =  b'  for  b  by  forward 

substitution. 

3  . 

Solve 

L’1BBTL'Ty  =  b  for  y 

by  conjugate 

gradient  iteration. 

4. 

Determine  x  by  back  substitution 

in  LTx  =  y. 

1 .  The  Incomplete  Factorization 


The  procedure  presented  here  is  given  in  Ajiz  and  Jennings 
[Ref.  13:  pp.  951-952],  and  Jennings  and  Malik  [Ref.  14:  pp. 
310-313].  To  see  how  the  elements  in  column  j  of  L  are  computed 
consider  the  following.  From  matrix  equation  3.12  a  typical  elemental 
equation  may  be  written  as 


1 . .  1 . .  r  p . ,  ♦  r» 

jj  >j  ij  ij 


]-i 


k..i 


likJjk 


( i— j ) 


(3.13) 


where  the  subscripts  for  elements  1  refer  to  their  positions  in  the  lower 

>T 

triangular  factor  L,  and  the  e jj  are  the  elements  of  BB  .  Let 


6  «  =  e«  'k?i  liklik  • 


(3.14) 


Then,  assuming  the  elements  of  L  have  been  computed  for  columns  1  to 

}|c 

(j-1),  all  the  e  in  column  j  can  be  computed.  First  consider  the  case 

a,c  • 

where  all  the  e  y  pass  the  rejection  test,  i.e.,  are  to  be  retained. 
Hence,  by  setting  Cy=0  in  equation  3.13,  we  get 


Xij  =  6  ' 


(3.15) 


j;c 

In  case  any  e  y  is  to  be  rejected,  ljj  is  set  to  zero,  implying  Cy=-e  jj  ■ 
The  off-diagonal  term  Cy  implies  diagonal  additions  Cy  and  c^-  to  the 
matrix  C  in  order  to  have  the  positive  semidefinitness  property.  The 
matrix  C  doesn't  have  to  be  stored  in  memory;  only  the  diagonal  addi¬ 
tions  Cy  are  of  further  interest. 

Before  any  ly  can  be  computed,  ly  has  to  be  determined.  With 
all  the  diagonal  additions  Cy  that  resulted  from  rejections  of  e'y  during 

the  computation  of  columns  1  to  (j-1),  an  expression  for  1,,  becomes 
-  JJ 


(3.16) 


where  e'y  is  defined  by 


p  =  p..  +  5* 

jj  jj  v  .  jj 


j-i 

2  1^  , 

(C*  l  1K 


(3.17) 


and  Cji^k^  is  the  diagonal  addition  to  C;;  resulting  from  deletion  of  e  ,  ■ , 

**  /  I  V  **  ...  J 

and  Cjd  '  the  diagonal  addition  to  c-.  resulting  from  deletion  of  e  ■  ^ . 

The  rejection  operation  tests  the  magnitude  of  an  element  e  y 
in  relation  to  the  current  values  of  the  corresponding  diagonal  elements 
ey  and  ejj  respectively,  whose  values  are  gi-ven  by 


where  the  index  k  refers  to  the  rows  for  which  rejections  in  column  j 

have  had  a  bearing  on  eH  and 

JJ  » 


e  ■  =  e-.  +  Y  c(k> 
11  11  tc-i  11 


(3.19) 


In  equation  3.19,  k  refers  to  the  columns  for  which  rejections  in  row  i 

_  >{c 

have  had  a  bearing  on  e^.  An  element  e  -  is  rejected  if 


*2  ,  2  -  - 

e  ij  <  *  ejjeii  * 


(3.20) 


where  ^  is  the  preset  rejection  parameter.  A  choice  of  ^=0  will  retain 
all  elements,  thus  leading  to  a  complete  Cholesky  factorization.  A  choice 
of  \p—  1  will  cause  all  off-diagonal  elements  to  be  rejected.  Ajiz  and 
Jennings  [Ref.  13:  p.  952]  recommend  that  the  rejection  parameter  be 
in  the  range  0.01<(p<0.2  for  effective  incomplete  factorizations. 

The  diagonal  modifications  in  C  which  result  from  the  rejection 
of  an  element  e '  y  are  given  by 


°jj(k)  =  ie\]i<ejj/ekk> 


(3.21) 


and  by 


i(k)  =  le  ikl  (eii/ekk) 1/2 


(3.22) 


With  the  successive  application  of  equations  3.14  in  conjunction  with  the 
rejection  operation,  3.17,  3.16  and  3.15,  all  elements  in  column  j  of  L 
are  determined. 

2 .  The  Conjugate  Gradient  Iteration 

The  conjugate  gradient  method  of  Hestenes  and  Stiefcl 
[Ref.  18]  is  applied  to  matrix  equation  3.10.  It  uses  the  following 
vectors,  the  letter  k  indicates  the  k-th  iteration, 

a)  p'  '  conjugate  gradient  vector 


b)  r' 


residual  vector 


c)  y(^)  solution  vector 

d)  u<k)  product  of  L"1BBTL'T  and  . 

The  initial  values  for  k=0  are  y^  =  0  and  =  r^)  =  b.  The  algo¬ 

rithm  for  one  iteration  is  as  follows 

u(k)  =  L-1BBTL“Tp(k) 
dk  =  (r(k))Tr(k)/(p(k))Tu(k) 

y(k+l)  =  y(k)  +  akp(5<) 

r(k+l)  =  r(k)  .  *  u(k) 

TV 

j3  =  (r(k+l))Tr(k+l)/(r(k))Tr(k) 

p(k+l)  =  r(k+l)  +  /3kP(k)  _ 

The  first  step  in  the  above  algorithm  is  obtained  without  computing  the 
transformed  matrix  explicitly  by  the  following  three  operations, 


1. 

LTv(k) 

=  p(k) 

(back  substitution) 

2. 

w(k> 

=  BBTV^k^ 

(pre -multiplication) 

3  . 

Lu(k> 

=  w  ^ k  ^ 

(forward  substitution) 

The  algorithm  is  terminated  when 

l|r(k)|l  /  llb||  <  tolerance  (3.23) 

where  b  equals  the  starting  residual  ,  since  y^  was  chosen  to  be 


zero. 


IV.  MODIFICATIONS  AND  COMPUTATIONAL  RESULTS 
A.  ALGORITHMIC  MODIFICATIONS 

1.  Iterative  Improvements 

A  provision  to  improve  the  solutions  to  equations  2.21  and 
2.22  has  been  implemented,  because  their  residuals  are  often  unduly 
high.  The  residuals  are  computed  as 

r  =  B1B1TA  *  b'  ,  (4.1) 

where  A  stands  for  Aq.  Aj  and  A2  respectively  and  b'  is  generic  for 
BlV,  (B2)i  and  (B2)2. 

If  |  r j |  >  e  for  some  i,  where  e  is  set  to,  say,  10  ’  the 

following  systems  are  solved  for  A' 

=  -r.  (4.2) 

An  improved  solution  is  then  obtained  by  adding  A  and  A' 

B1B1t(A+  A')  =  r+b'+c-r)  =  b'.  (4.3) 

The  improvement  can  be  repeated  if  the  residuals  of  equation  4.3  are 
still  found  to  be  too  high.  With  this  modification  the  direct  method 
(MDOC)  has  become  a  semi- iterative  method. 

2.  Removal  of  the  Artificial  Column 

The  update  of  the  artificial  column  with  the  residuals  of  the 
current  solution 


A(x',x'n  +  1) 

‘  bx'n  +  2  =  r 

(4.4) 

-  b  = 

AX(k*1)  -  b  -  r/x'n  +  2(k) 

(4.3) 

35 


on  each  iteration  has  been  augmented,  such  that  the  artificial  variable 
can  be  deleted  from  the  problem  when  the  infeasibility  becomes  small. 
After  the  deletion,  only  one  dense  column  remains  in  B,  and  the  solu¬ 
tion  of  the  least- squares  problem  is  simplified.  The  motivation  for  this 
modification  is  to  avoid  some  of  the  numerical  problems  that  arise  from 
variables  that  approach  zero. 

3.  Positive  Semidefinite  BB? 

It  can  be  shown  that  if,  and  only  if  a  diagonal  element  ever 

T 

goes  to  zero  in  a  Cholesky  factorization,  BB  is  positive  semidefinite, 
not  positive  definite.  Furthermore,  all  elements  below  the  zero  diagonal 
element  must  also  be  zero  and  a  nonunique  solution  to  the  normal  equa¬ 
tions  can  be  obtained  by  setting  the  corresponding  A-  to  0.  This  solu¬ 
tion  can  be  obtained  by  replacing  the  zero  diagonal  element  with  any 
positive  value  and  continuing  with  the  factorization.  Restated,  if  a 
diagonal  element  1--  in  the  partially  computed  Cholesky  factor  is  less 

than  or  equal  to  10  ,  set  ljj= 1 ,  and  l-j=0  for  i=j  +  l, j  +  2,  .  .  .  ,m  in  equa¬ 
tions  3.2  and  3.3  respectively.  This  procedure  is  also  useful  to  deal 

with  numerical  problems  which  arise  from  degeneracy  near  optimality. 

4.  Weighted  Homogeneity  Variable 

The  initial  solution  (y° , y°n  +  1 , y°n+2)  =  1  used  to  begin  the 
projective  algorithm  is  arbitrary,  and  in  some  sense,  the  "homogeneity 
variable"  y°n  +  2  is  fundamentally  different  than  the  other  variables. 
Consequently,  a  modification  has  been  made  to  allow  weighting  the 

starting  solution  y°n  +  2  differently  from  the  other  y°.,  i=l,...,n  +  l.  Let 
y°n  +  2=s-  With  ir^'y°=n+2,  the  other  y°j  are  set  equal, 

y°i  =  (n+2-s)/(n+l) ,  i=l , . . . , n+1.  (4.6) 

It  is  hoped  that  such  a  weighting  might  lead  to  lower  iteration  counts. 
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B.  TEST  PROBLEMS  AND  COMPUTATIONAL  RESULTS 


At  first,  a  small  test  problem  (TEST1)  was  used  to  verify  the 
correctness  of  the  source  code.  The  two  algorithms  were  then  tested  on 
seven  problems  which  have  been  also  used  for  testing  by  Lustig 
[Ref.  4:  p.  19].  Table  5  shows  some  of  the  characteristics  of  the  test 
problems . 


TABLE  5 

TEST 

PROBLEMS 

Name 

Rows 

Columns  Logicals 

Nonzeros 

Density 

TEST1 

6 

2 

6 

10 

0.833 

AFIRO 

27 

32 

19 

95 

0.103 

ADLITTLE 

56 

97 

41 

522 

0.096 

1 

SHARE2B 

96 

79 

83 

901 

0.119  ! 

ISRAEL 

174 

142 

174 

2529 

0.102  1 

BRANDY 

193 

249 

54 

2204 

0.046 

.  E226 

223 

282 

190 

2578  . 

0.041 

BANDM 

305 

472 

0 

2.659 

0.018 

Table  6  summarizes  the  computational  results.  All  CPU  times 
represent  the  time  in  seconds  to  set  up  the  •  major  part  of  the  data 
structure,  solve  the  LP,  and  write  out  a  few  parameters  on  each  itera¬ 
tion  and  the  solution.  Times  to  read  in  the  data  from  MPS  format  are 
not  included. 

The  convergence  criterion  equation  2.31  is  used  with  t=0.05  for  all 
test  problems  but  TEST1  and  SHARE2B,  where  t=0.001  and  t=0.01 
respectively.  No  solutions  have  been  obtained  for  the  problems  A  TIRO, 
BRANDY  and  BANDM:  the  algorithm  will  not  converge  to  the  optimum. 

The  feasibility  parameter  p  (equation  2.7)  is  set  to  0.9995,  giving 
the  best  overall  performance  of  the  algorithm.  Tests  with  p=0.9999  and 


p=0. 99999  on  the  data  sets  TEST1  and  SHARE2B,  indicate  that  with 
these  higher  values  the  number  of  iterations  necessary  to  attain  feasi¬ 
bility  is  reduced,  but  the  total  number  of  iterations  remains  about  the 
same.  In  addition,  the  accuracy  of  the  solutions  deteriorates. 

An  intermediate  solution  is  declared  feasible,  i.e.  the  artificial 
variable  is  removed  from  the  problem,  as  soon  as  xn+jM  is  less  than 
10  ^ .  The  choice  of  10  ^  is  arbitrary,  and  depends  on  the  value  of  the 

_  O 

optimal  solution.  Using  alternate  values  of  1.0  or  10  makes  little 
difference  in  the  number  of  iterations  necessary  to  attain  feasibility. 


TABLE  6 

COMPUTATIONAL  RESULTS 
MDOC  ICCG 


Problem 

Iterations 

feasible  optimal 

Time 

Iterations 

feasible  optimal 

Time 

TEST1 

1 

6 

0.03 

1 

6 

0.04 

ADLITTLE 

13 

21 

1.96 

13 

21  . 

2.18 

SHARE2B 

6 

27 

5.04 

7 

37 

5.31 

ISRAEL 

19 

98 

311 

19 

89 

347 

E226 

8 

49 

107 

9 

49 

130 

Factorization  failures  plagued  both  algorithms  before  the  modifica- 
tion  to  accomodate  a  semidefinite  BB  was  implemented.  These  failures 
occured  near  the  optimum,  e.g.  E226,  or  with  the  ICCG  algorithm  when 
setting  the  rejection  parameter  to  a  value  greater  than  zero.  After 
implementing  the  semidefinite  modification  these  factorization  problems 
have  been  cured,  but  the  ICCG  algorithm  now  shows  very  slow  conver¬ 
gence.  For  example,  a  solution  to  SHARE2B  is  obtained  only  after 
112.36  CPU  seconds  with  \p- 0.015  and  all  other  parameters  unchanged. 
Thus,  because  storage  is  not  at  a  premium,  only  complete  factorizations 


are  used.  Satisfactory  numerical  results  with  the  incomplete  factoriza¬ 
tion  are  obtained  only  with  the  trivial  test  problem  TEST1;  computation 
times  are  always  inferior. 

The  minimum-degree  ordering  works  very  well  in  practice. 

Computational  cost  seems  to  be  moderate,  e.g.,  0.25  seconds  for 

SHARE2B,  and  3.12  seconds  for  BANDM.  Table  7  gives  densities  of 
T 

BB  and  the  Cholesky  factors  with  and  without  reordering.  Substantial 
reductions  in  storage  requirements  were  achieved  when  doing  the 
reordering. 

TABLE  7 
DENSITIES 


Re-ordered  Not  re-ordered 

Nonzeros  Nonzeros  Nonzeros 


Nome 

in  BB 1 

Density 

in  L 

Density 

in  L 

Density 

TEST1 

20 

0.952 

20 

0.952 

20 

0.952 

AFIRO 

62 

0.177 

107 

0.305 

194 

0.553 

ADLITTLE 

384 

0.241 

411 

0.258 

816 

0.512 

SHARE2B 

871 

0.187 

1021 

0.219 

1134 

0.243 

ISRAEL 

11227 

0.737 

11433 

0.751 

13743 

0.903 

BRANDY 

2853 

0.152 

3429 

0.183 

9760 

0.521 

E226 

2823 

0.113 

3639 

0.146 

10735 

0.430 

BANDM 

3724 

0.080 

4660 

0. 100 

32090 

0.688 

The  densities  and  number  of  nonzeros  in  L  or  BB  are  relative 
to  a  symmetric  half  of  a  matrix;  the  diagonal  elements  are  not 
included . 


choice  of  M.  Doing  a  few  (2  or  3)  iterative  improvements  at  this  stage 
stabilizes  the  computations  until  the  algorithm  gets  close  to  the  optimal 
solution.  Numerical  problems  near  the  optimum  are  so  severe  that  even 
allowing  a  prohibitively  high  number  like  25  iterative  improvements  has 
no  apparent  influence  on  the  quality  of  the  solution. 

The  mildest  form  of  numerical  difficulty  near  the  optimum  is  slow 
convergence,  e.g.  ISRAEL  and  E226.  For  SHARE2B,  a  high  quality 
solution  is  obtained  with  no  numerical  difficulty.  This  data  set  is  chosen 
to  test  the  weighted  homogeneity  variable  modification.  S  is  given 
several  values  in  the  range  1  to  100.  Iteration  counts  range  from  25  to 
37.  The  quality  of  the  solutions  is  sometimes  reduced,  however.  Only 
s=35  (27  iterations)  and  s=50  (28  iterations)  give  high  quality  solutions 
with  low  iteration  counts.  Values  of  s  greater  than  75  creat  conver¬ 
gence  failures. 

C.  CONCLUSIONS 

The  low  iteration  counts  for  some  of  the  test  problems  are  prom¬ 
ising,  although  CPU  times  seem  to  tell  the  difference.  These  high  CPU 
times  result  from  test  problems  having  slow  convergence  near  the 
optimum,  which  is  believed  to  be  due  to  many  variables  going  to  zero. 
Thus,  a  technique  to  drop  variables  going  to  zero  from  the  LP  could 
well  speed  up  convergence. 

The  ICCG  algorithm  does  not  perform  very  well  on  the  test  prob¬ 
lems.  Computational  results  are  generally  inferior  to  those  obtained 
with  the  MDOC  algorithm.  Thus,  no  further  research  into  this  method 
seems  warranted. 

In  view  of  the  numerical  problems  encountered  when  solving  the 
least-squares  problem  with  the  normal  equations  approach  and  Cholesky 
factorization,  another  method  that  does  not  use  square  roots  should  be 
considered  for  implementation.  A  very  promising  candidate  in  this 
respect  is  Givens  rotation.  See  Gentleman  [Ref.  19]  and  George  and 
Heath  [Ref.  20]. 
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